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Abstract 


We calculate the maximal Lyapunov exponent in constant-energy molec¬ 
ular dynamics simulations at the melting transition for finite clusters of 6 to 
13 particles (model rare-gas and metallic systems) as well as for bulk rare- 
gas solid. For clusters, the Lyapunov exponent generally varies linearly with 
the total energy, but the slope changes sharply at the melting transition. In 
the bulk system, melting corresponds to a jump in the Lyapunov exponent, 
and this corresponds to a singularity in the variance of the curvature of the 
potential energy surface. In these systems there are two mechanisms of chaos 
- local instability and parametric instability. We calculate the contribution 
of the parametric instability towards the chaoticity of these systems using a 
recently proposed formalism. The contribution of parametric instability is 
a continuous function of energy in small clusters but not in the bulk where 
the melting corresponds to a decrease in this quantity. This implies that the 
melting in small clusters does not lead to enhanced local instability. 

I. INTRODUCTION 

In recent years the dynamics of finite condensed matter systems, especially atomic and 
molecular clusters, has been extensively studied from a nonlinear dynamics perspective [|I|. 
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Various quantities of interest like Lyapunov exponents, Lyapunov spectra, distributions of 
finite-time Lyapunov exponents and the Kolmogorov entropy have been computed to see 
evolution of chaoticity and ergodicity. Very general considerations lead to the expectation 
that the Lyapunov exponents and the Kolmogorov entropy should increase with energy. 
However, there are indications that different systems can display a variety of behaviour, and 
details of how such indices change and the different kinds of possible (qualitative as well as 
quantitative) behaviour—the various universality classes, so to speak—have not yet been 
completely characterized. 

In the present work we calculate the largest Lyapunov exponent. A,for small rare-gas 
and metal clusters (modelled, respectively, by the Lennard-Jones (LJ) and the many-body 
Gupta (Gu) 0 potentials) as well as for bulk rare-gas solid, namely 256 LJ atoms in a box 
with periodic boundary conditions. In the range studied, A is generally linearly related to 
energy (except at very low energies) but shows a sharp change in slope at an energy which 
can be related to the melting transition. 

We also compute an estimate for A through a semi-empirical methodology provided by a 
recently proposed geometrical theory of hamiltonian chaos Q. Under certain approximations 
this allows for the estimation of the relative contribution of stable modes of a hamiltonian 
system towards chaotic behaviour. The approximations inherent in the theory P] are fulfilled 
in bulk systems (where N is large), but do not appear entirely satisfactory in small clusters. 
This geometrical theory has as its input the curvature of the configuration space manifold 
and its fluctuation. These quantities and their variation with temperature are themselves 
interesting because they yield a statistical quantification of the potential energy surface. 
Recent approaches to the computation of Lyapunov exponents from (local) instantaneous 
mode analysis use this information implicitly. 

The simplest understanding of chaotic dynamics in such (hamiltonian) systems is the 
standard KAM picture |p. When the number of freedoms becomes large (for 13 atoms, 
the phase space is 78 dimensional) the picture of a phase space foliated by tori, with sur¬ 
rounding stochastic regions p is not particularly relevant. However, the motion for specific 
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initial conditions remains nearly periodic, while for others there can be a positive Lyapunov 
exponent. In particular the KAM theorem underestimates the chaotic thresholds by several 
orders of magnitudes |]^ in high-dimensional systems as the critical perturbation scales as 
~ exp —aNf which rapidly goes to zero with increasing Nf (degrees of freedom), contrary 
to the experience in numerical simulations beginning with FPU’s famous result [^. And 
attainment of chaoticity does not exhaust the interest in dynamics—in particular the evolu¬ 
tion of the dynamics near a thermodynamic phase-transition is nontrivial. At the energies 
corresponding to the phase-transition phenomenon, the accessible phase space increases 
greatly, and correspondingly A shows a signature of the transition. An alternate means of 
analysis is through decorrelation of the eigenvectors of the instantaneous Jacobian along a 
trajectory Q. The timescale for this process is greatly reduced by the presence of the regions 
of negative curvature (unstable modes) which also increase at the melting phase-transition 
point. 

Such ideas have been at the root of a variety of studies of the Lyapunov exponent or 
related quantities. Posch and Hoover calculated Lyapunov spectra of solid and liquid 
LJ systems in two and three dimensions, and attempted to £t a power law to this data, 
obtaining different exponents in the solid and liquid phases, although no dehnite significance 
could be ascribed to these. Berry and coworkers m have examined a variety of quantities 
including hnite-time or local approximations to Kolmogorov entropy and Lyapunov spectrum 


in Lennard-Jones (LJ) and Morse clusters containing between 3 and 13 particles These 
studies have been able to make contact between the features of the potential energy surface 
and the variation of different dynamical indicators. Recently Sastry ^ has computed the 
maximal Lyapunov exponent as well as the entire Lyapunov spectra for a 32 atom Lennard- 
Jones fluid in the temperature range 50-800 K. A power law £t of A with temperature does 
not yield a single exponent but a crossover (around T = 1) between the exponent = 1 at low 
temperature to 1/2 at high temperatures |]^. In a hnite lattice system, Butera and Caravati 
0 simulated the coupled rotor system in two dimensions which has a Kosterlitz-Thouless 
(KT) transition at constant energy and observed that the scaling of A with the temperature 
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changes at precisely the KT temperature. Recent simulations on the XY model in 2 and 
3 dimensions by Caiani et. al. 0 have shown the differences between KT transition and 
true symmetry-breaking transition in three dimensions. The crossover in scaling of A with 
temperature or energy per particle, also observed in other lattice models where it had purely 
dynamical significance |^, was suggested to coincide with the crossover from slow to fast 
diffusion in the phase space [0 and was labelled as the strong stochasticity threshold (SST). 
The generality of SST in nonlinear hamiltonian systems is not obvious, although it appears 
to persist even for large degrees of freedom in the models studied. It was also found that 
SST occurs in lattice models with Lennard-Jones interactions in two and three dimensions 
, the signature being a rapid relaxation of the specific heat and independence of A on the 


initial conditions. Similar transitions have also been observed in small rare-gas and metal 


clusters 0]. A somewhat different interpretation of these results has also been proposed 

ra¬ 
in small rare-gas atomic clusters, A has been calculated at the melting transition ||19|| , 
which is a hnite-size analogue of the bulk melting transition [0|. Whereas A changes 
discontinuously with energy for LJ 13 and LJ 55 , for LJy, the slope changes discontinuously. 
At the energy of this change, indicators of melting like Lindemann index or density of states 
show characteristic signature, so it was proposed that A is a good indicator of melting 
transition 0,0]. Subsequent work p 2 [| has revealed that the behaviour of A as a function 
of the energy is more complicated and can be different depending on the nature of the low 


energy configuration that the system starts from. More recently Dellago and Posch |23 


have calculated A,the Lyapunov spectra and related quantities like the fraction of unstable 
modes for the melting transition of modified LJ systems in 2D (the so-called WCA and LJS 
potentials). They hnd that A has a broad peak near the melting density and a steplike 
increase at the melting temperature, and further suggest that there is a change in the 
shape of the Lyapunov spectra at the transition density. Critical phenomena which occur 
at higher temperatures in larger fragmenting clusters have also been studied 0 ], although 
the claimed form for the (finite-time) Lyapunov exponent, namely A ~ (T — Tc)~^, with 
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universal behaviour for /i is questionable HI^,^ . 

Our work in the present paper is focussed on the study of clusters and bulk at the 
melting transition, with particular emphasis on the Lyapunov exponent and the curvature 
fluctuations. In the next section, we examine the behaviour of A,and in Section 3, the 
nature of the curvature fluctuations are analysed in order to make contact between theory 
and simulations. The methodology and theoretical background for the extraction of A from 
curvature data is elaborated in Section 3, where our results for bulk as well as cluster systems 
are also presented. This is followed by a summary and discussion in Section 4. 


II. MELTING AND LYAPUNOV EXPONENT IN CLUSTERS AND BULK 


As is by now well-known atomic clusters with as few as six or seven particles 
undergo a hnite size analogue of the bulk melting phase transition, which is marked by a 
jump in the Lindemann index and the onset of rapid isomerization. The liquidlike phase 
of the small clusters is however unlike the bulk liquid in one important respect: due to 
the phase space constraints, the atomic displacements (“hopping processes”) are highly 
correlated, giving rise to 1// spectra of single particle potential energy fluctuations |^. It 
has been observed that in metal clusters particles can occupy two types of sites—of low and 
high energy respectively, and the onset of the liquid phase corresponds to the onset of an 
isomerization occuring through four particle interchange between high and low energy sites 
Similar features are expected for rare-gas clusters also. This peculiar dynamics has 


interesting consequences which we examine in this section. 

We consider clusters of up to 13 atoms interacting via Lennard-Jones (LJ) potential 
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which is appropriate for rare-gas atoms; we work in reduced units when cr = e = 1. In order 
to model metallic clusters the manybody Gupta (Gu) potential |^, 

V = (l/2)t/^{A^exp(-p(rij -ro)) 
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-[^exp(-2g(ry -ro))]^/^}. 


(2) 


is commonly used. Here rij is the distance between the atoms i and j and ro is the interatomic 
distance in the bulk (fee) crystal. The specihe metallic cluster system we model corresponds 
notionally to Ni, for which we use parameters given in Ref. [Q: p = 9/rQ,q = 3/rQ,A = 
0.101035 and U/E^uik = 0.85505 {Ebuik is the bulk cohesive energy of the metal by which all 
energies are scaled). Quantities like A and Kolmogorov-Sinai entropy have been calculated 
for small rare-gas clusters and the relation of the potential energy surface to local dynamical 
behaviour has been analysed ||l0|,[l^. In particular it is known that smooth saddles cause a 
drop in local chaoticity indicators . 


Time is measured in units of (ma^/e)^/^ for LJ clusters and {mrl/EbuikY^^ for Gupta 
clusters. Constant energy simulations are done using the Verlet algorithm with stepsize 
At = 0.01 in these units, and the total energy is conserved to 0.01%. The total integration 
time varies between 10^At to 5 x 10®At depending on the potential and the system size. 
Simulations start at the global minimum of the structure and then gradually move up in 
energy; at the highest temperatures studied evaporation does not set in within the time 
period of the computation. Temperature is dehned in the usual manner, as proportional to 
the average kinetic energy per particle, T = 2{Ek)/{kB{3N — 6 )), ks being the Boltzmann 
constant. 

At very low energies the dynamics is nonergodic, in particular for 13 particle clusters the 
nonergodicity can be persistent upto rather high energies. The breathing mode in LJ 13 was 
recently studied by Salian et. al. |^. This mode is stable (h e. nondecaying) up to an 
energy E = 0.150 per particle above the global minimum; the corresponding mode for LJ 7 
survives only up to energy 0.042 above the minimum. Stability here is tested by starting 
trajectories with isotropic stretching of the global minimum structure. In this nonergodic 
region the system initialised with small random distortion of the icosahedral structure has 
a very small positive Lyapunov exponent while for larger distortions at the same energy, 
A can be much larger. As energy increases, the overall cluster distortions increase and this 
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mode becomes markedly unstable. For large initial distortions from icosahedral geometry 
of LJi 3 , a chaotic transition occurs at a lower energy and is manifested as the divergence 
of the microcanonical specihc heat. This energy depends on the equilibration time but 
tends to a nonzero value for large equilibration times. Such transitions are not size-specihc, 
and have been seen for non-magic clusters as well. 

Fig. 1 shows the variation of A(e) with the reduced energy per particle, e, for LJ and Gu 
clusters of various sizes. At higher energies A(e) is linear but at a certain energy a sharp 
change in the slope is evident in all cases. Precisely at this energy the conventional criterion 
of melting applies, namely the Lindemann index crosses the value 0.1 (Lindemann indices 
for Gu clusters have been studied in |Q). In LJ 13 which has a large solid-liquid coexistence 
region between —2.6 < e < —2.2, A changes slope at e ~ —2.6. In LJg where the Lindemann 
index increases continuously, the discontinuity in A appears just after the Lindemann index 
has reached its critical value 0.1. The slope in the liquid phase is always smaller, indicating 
that the diffusive modes have different chaotic properties, giving rise to different energy 
dependence of A.The change in the slope of A(e) can also be taken to be a characteristic 
signature of melting in small clusters. The slope of the A(e) curves are generally smaller for 
the larger cluster; furthermore, the sharpness of the discontinuity (in slope) is reduced as 
the cluster size increases. 

For the Gupta clusters, however, these trends with N are not strongly-marked which is 
a consequence of the manybody character of the Gu potential: even if a pair of particles is 
not interacting directly (being outside the potential cutoff) the corresponding elements of 
the Hessian matrix need not vanish because of the manybody term. We hnd that the slope 
changes distinctly at an energy corresponding to top of the jump in the Lindemann curve. 

The third system we study is bulk, and Fig. 2(a) shows A(e) for the system of 256 LJ 
particles in a cubic box with periodic boundary conditions at the reduced density p = 0.93 
and reduced melting temperature 1.15. Initial conditions for these simulations were as 
follows: the atoms were initially at fee lattice positions, with initial velocities taken randomly 
from an appropriate gaussian distribution. For e > —4.25 the crystal is unstable and soon 
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melts. It is possible for melting to occur for slightly lower energies if integration is carried 
for longer times but the time required to melt is not a monotonicahy decreasing function 
of energy. Here A(e) shows a jump which obviously can be ascribed to the increase in the 
disorder at melting. The data shown is the Lyapunov exponent of solid phase or the liquid 
phase trajectory only which are obtained by discarding the premelting portion of a trajectory 
that melts. Similar results for bulk melting have also been reported by Dellago et. al. 


and Nayak et. al. 


While both the bulk and the cluster are disordered by melting (change in entropy at 
melting, AS/N = 1.0 for LJ 55 and 1.7 for LJ crystal and specihc heat of even 6 particle 
cluster shows a peak) A in the cluster liquid phase is signihcantly smaller than the value 
obtained by extrapolating from the solid phase, in contrast with A obtained in bulk liquid. 
This suggests that in the cluster liquid phase there are stabilizing influences on the dynamics 
which are absent in the bulk liquid. We conjecture that the correlated hopping processes 
in clusters [^,^] provide the necessary mechanism. This is consistent with the observation 
made above that the A(e) curve gets smoother with the cluster size. 

The observed Lyapunov phenomenon for the clusters is not just an effect of the smeared- 
out bulk transition i. e. a manifestation of hnite dynamical coexistence region in the clusters 
which vanishes in the bulk limit. The properties of the coexistence region (in particular its 
width) depend sensitively on the cluster geometry e. g. the fact that the cluster is magic 
or not. But the trends in the Lyapunov exponent depend on the size in a simple manner 
and are independent of the magic-nonmagic relation. 

The spectra of 3N — 7 positive Lyapunov exponents for the clusters vary smoothly with 
energy. For LJ clusters they are linear in the entire range with a slope that increases with 
energy while in Gu clusters they acquire increasing curvature. This is contrast to the results 
of Dellago and Posch for bulk melting in two dimensions where the curvature of the 
spectra changes sign smoothly at certain densities (at constant temperature). It remains a 
task to extract more useful information from the shape of the Lyapunov spectra. However 
the smoothness of the spectra at cluster melting implies that the relation (if any) between 
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thermodynamic and dynamic entropies is nontrivial. 


III. CURVATURE FLUCTUATIONS 


In this section we apply the geometric theory of Hamiltonian chaos developed by Pettini 
and coworkers p] in order to interpret and rationalize the results of our numerical simulations 
in terms of an underlying (microscopic) mechanism. This theory, the salient features of which 
are summarized below, is attractive because it attempts to unite features of the potential 
energy surface with the dynamical properties of the system as encoded in the Lyapunov 
exponents. One additional motivation in applying this theory to hnite cluster systems is 
to determine the limits of applicability of the general framework, which has mainly been 
applied to lattice models where the calculated Lyapunov exponents are found to be in good 
agreement with empirical exponents P,P . 


A. The Geometric theory of chaos in High-dimensional systems 


It is well known that the classical dynamics on manifolds of constant negative curvature 
is chaotic . The dynamics on a manifold with fluctuating positive curvature can also be 


chaotic this fluctuation, via the mechanism of parametric instability, is responsible for 
creating chaos in systems such as the Fermi-Pasta-Ulam (3 and 0"^ chains 0 and for a driven 
one-dimensional oscillator studied by Chaudhuri et. al. [0]. These studies have provided 
much of the motivation for the development of the geometric theory of chaos p. Barnett 


et. al. have applied similar ideas to calculation of Lyapunov exponent of a dilute gas. 

The geometric theory makes a diagonal approximation in the sense that it uses infor¬ 
mation only about the trace of the instantaneous Hessian matrix i. e. AV. When 
the equations of motion are put in a differential-geometric form, this term appears as the 
Ricci curvature (curvature locally averaged over the directions) of the enlarged conhguration 
manifold in the Eisenhart metric |0|, 
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ds'^ = —2V {q)dt^ + aijdq^dq^ + dtdq^^^^. 


( 3 ) 


Here V{q) is the potential energy and the kinetic energy is T = atjq^qK The relevance 
of potential curvature has been noted before [^. The Ricci curvature does not have this 
simple form in other metrics but the essence of the assumption is that all the directions in the 
phase space are equally curved after a coarse-graining along a trajectory. The dynamical 
trajectories are the geodesics of this manifold. In above models this appears to be the 
dominant mechanism for chaos as there are no unstable modes (corresponding to the negative 
eigenvalues of the instantaneous Hessian matrix or regions of negative curvature) which are 
local mechanism of chaos. 

If it is assumed that the curvature fluctuations have gaussian spectra upto a high- 
frequency cutoff i. e. the dynamics generates a gaussian random process for curvatures, then 
one can dispense with the necessity of following a trajectory and an estimate of the largest 
Lyapunov exponent, A can be obtained via Monte Carlo sampling. This is the essence of the 
gaussian approximation 0, within which excellent agreement is obtained between A and A. 
The presence of additional unstable modes renormalizes the calculated exponent, although 
this is nontrivial to calculate. Therefore the nnr enormalized exponent gives an estimate of 
the chaoticity caused by stable modes only (the unstable modes also contribute to paramet¬ 
ric instability by their presence in AV but it is not their major effect on chaoticity). The 
theoretical basis of the diagonal approximation assumption is that the local fluctuation of 
the Ricci curvature detects deviation from isotropy (at a point) because isotropic manifolds 
are necessarily of constant curvature by Schur’s theorem |P]. 

The crossover between the regimes of weak and strong chaos in high-dimensional systems 
can be detected by examining the behaviour of the mean curvature, A;, as a function of the 
energy density. In the integrable limit k is independent of energy [^. Corresponding to a 
Hamiltonian H oi N particles with an interaction V, 




( 4 ) 


there are 3N equations of motion. 
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The associated Jacobi equation for the second deviations is then 




dt"^ 


dqidqj 


( 6 ) 


where J* are the components of the vector of second deviations. After some approximations 
PJ this can be converted to an equation for u = \J\, 


d?u d‘^V JiJj d?u , . 


( 7 ) 


which is, in effect, an equation for a linear oscillator with time-dependent frequency \fQ. 
The solutions of this equation are unbounded for typical Q{t) and the Lyapunov exponent is 
just given by the rate of exponential growth of the envelope of u Pettini et. al. P show 
that Q{t) is just the sectional curvature (which is the generalization of gaussian curvature 
to many dimensions) relative to the plane containing J and dq/dt in the Eisenhart metric. 
The diagonal approximation consists in replacing Q by the simpler quantity 


AE _ 1 _ d‘^V 


( 8 ) 


with Nf the number of degrees of freedom, which is p, the Ricci curvature per degree of 
freedom i. e. sectional curvature has been averaged over relative orientations of J and the 
velocity vector. Under the further assumption that the curvature is gaussian distributed 
with a mean k = {AV)/Nf and variance = (((AU)^ — {AV)‘^)/Nf and are J-correlated, 
Pettini et. al. p derive an expression for an estimate of the Lyapunov exponent 


A = i(/ 




l = a^T + 


IQAk^ 


1/3 


_|_ ((j27-j2 


( 9 ) 

( 10 ) 


where r is a characteristic time implied by the smoothness of the underlying manifold i. 
e. the time-interval below which dynamics of curvatures cannot be regarded as a random 
process. 
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B. Application of the Geometric theory 


The main result of the geometric theory is an estimate of the Lyapunov exponent, A 
given in Eq. (|^), for which it is necessary to obtain the mean curvature, k and the variance, 
cr^. These quantities can be calculated along a typical trajectory using Eq. (H), and the 
assumption of 5—correlated curvature fluctuations can be directly verified. 

The determination of the timescale r (see Eq. 10) is more tricky. One estimate which 
has been used |^, 

Ti^/k 

2^ k(k + a) + TTcr 

This expression for r here is actually that for r* in Ref. [Q (see the discussion following 
Eq. 45 there). However in the presence of negative curvatures it may be more accurate to 
use a different timescale 0 

^2 = - 

a 

We find (see the next subsection) that T 2 is more accurate than r insofar as it provides a 
better numerical match with the autocorrelation decay timescale for the systems studied 
here. 

One additional minor point is that the effective number of freedoms is Nj = 3N — 6 
rather than 3N, since the six conserved quantities (linear and angular momenta) give rise 
to zero frequency modes and thus do not contribute to the chaoticity. This does not change 
qualitative conclusions (indeed it must not) and improves results in most cases. 

The application of the geometric theory to the systems considered in Sec.^ is of interest 
for two reasons. Firstly, this formalism has so far been mostly applied to lattice models, 
where parametric instability is the main source of chaos. It would be useful to determine the 
extent to which the formalism works for off-lattice models with significant local instability. 
Secondly, the partition of the chaoticity of the system into local and nonlocal components 
may help in clarifying the behaviour of Lyapunov exponent at melting. In particular, it 
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would be interesting to know whether the overall instability of the system can be separated 
into these two components, and if so, whether they behave differently at phase transitions. 

The melting transition in hnite systems appears to have a distinct signature—the Lya¬ 
punov exponent shows a knee, where the slope changes discontinuously [^. For bulk, the 
fraction of the unstable modes has a discontinuous jump at melting. Therefore a rapid 
increase in local instability and consequently, a jump in the Lyapunov exponent can be ex¬ 
pected. Such a change may not be apparent in the contribution of the parametric instability, 
and therefore in the situation of cluster melting where Lyapunov exponent does not increase, 
it is an open question whether the local instability increases or not. 


C. Results 

As should be clear from the preceeding discussion, application of the geometric theory in 
order to make comparison with the results of our numerical simulations involves a number 
of sensitive estimates and several approximations. 

Following the general procedure outlined in Sec. IIIB above, we have calculated the 
estimate A (cf. Eq. ||) for bulk (LJ) and various LJ and Gu clusters in an energy range 
which encompasses melting, from long trajectories of duration up to 2 x 10®At. The data 
for k, a and A is shown in Figs. 2-6. The detailed comparisons of theory and simulations are 
presented separately for bulk and cluster systems below. 


i. Bulk LJ system 


Casetti and Macchi have calculated the curvature for bulk LJ in a exponentially 
large energy range in order to detect the crossover between weakly and strongly chaotic 
regimes. Earlier calculations by LaViolette and Stillinger of the mean curvature (which 
is proportional to the squared Einstein frequency) show that k increases linearly in the 
bulk LJ system with a jump of about 20% at melting. This increase is a manifestation of 
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the positive high-frequency tail in the instantaneous normal spectrum in the liquid phase. 
However, the slope in the liquid phase is smaller than the solid phase by about 6%. 


We hnd that the variance cx^ has a discontinuity at an energy = —4.17. At energies 
well away from e^, cr increases linearly with energy but in a very narrow range preceding 
Cm, roughly corresponding to the solid-liquid coexistence region, a increases sharply. As 
the system melts a falls by about 30%. (The discontinuity in a has been confirmed by 
repeating the calculations with longer trajectories and hner energy mesh. Data in the 
coexistence region was computed from long trajectories of total time 2 x lO^At. Care was 
taken to ensure that computed averages are over either solid or liquid phase exclusively.) 
This behaviour may be contrasted with the cusp singularity found recently in XY model in 


three dimensions by Caiani et. al. which was interpreted as suggestive of a topological 
transition in the potential energy surface. 

The assumption of h—correlated curvature fluctuation can be substantiated by examining 
the power spectra of curvature fluctuations. Fig. 7 shows that these are satished in the bulk. 
However the observed correlation time does not agree with r given in Eq. o (see Table 1). 
We therefore use T 2 to calculate the relative contribution of unstable modes, namely 


<5A„ = (A-A)/A. (13) 

Indeed A(r 2 ) is much better £t to A than A(r) (Fig 2). 6Xu can become slightly negative in 
the solid phase (implying some correction due to correlations), while in the liquid phase, it 
can be as large as 0.35. Assignment of the difference A — A to the effect of unstable modes 
is justified by the following two observations. Firstly, A does not increase at melting (it 
actually falls) whereas A has signihcant jump which can accounted for by an increase in 
the fraction of unstable modes.. In addition, the agreement of A and A is better at lower 
temperatures, namely when the occurrence of negative curvatures is infrequent. 
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2. Clusters 


Owing to the finiteness of cluster systems, correlations do not decay sufficiently rapidly 
2 ^. As a consequence, curvature fluctuations are far from being uncorrelated, and the 


framework of the geometric theory breaks down for such systems and deviations from Eq. 
can be expected. The observed correlation times do not agree with analytical estimates for 
r and T 2 (Table 1). 

The mean curvature for specific LJ clusters has been computed previously 00 , and in 
contrast to bulk, decreases uniformly with energy for all clusters (except LJ 13 ); likewise for 
Gupta clusters. No trend is apparent either for k or its variation with energy, although there 
are some size effects in the case of Gupta clusters. The behaviour of the variance, is more 
complex. Although this quantity usually increases smoothly with energy, in the coexistence 
regime near melting there are large fluctuations which persist for very long averaging times. 
The liquid phase in LJ 13 also shows nonmonotonic dependence of variance on energy. 

The net result of the persistence of correlations is that the estimates A for cluster systems 
do not agree with A.While A(r) is a rather good fit for LJ clusters, it is of doubtful validity 
because r is far from the observed fluctuation timescale Tc (Table 1). Using T 2 in Eq. 
gives A/A 1.3 — 2.5 for LJ clusters, although for the tightly bound Gupta clusters, this 
discrepancy is smaller, \/Kk. 0.8 — 1 . 2 . 

3. Discussion 


Our results indicate that unstable modes have suppressed chaoticity in certain circum¬ 
stances. In particular, in the solid LJ system, where 6Xu is very small, the fraction of unstable 
modes is substantial (~ 0.2) while slightly higher (~ 0.25 — 0.3) fraction of unstable modes 
in liquid gives greatly enhanced 6Xu (~ 0.35). In the cluster even after melting 6Xu does not 
increase very much: it has smooth dependence on energy. A tentative conclusion that can 
be drawn from these cases is that the unstable modes have effective chaoticity only when 
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the particles are free to execute large-scale motion. 

If the coarse-grained curvature is everywhere positive, the ratio a/k provides a crude 
measure of the ruggedness or roughness of the underlying potential energy surface. As it 
is this feature which causes nearby trajectories to diverge, it is interesting to study the 
variation of this index with energy, as this will give some indication of the nature of the 
region that is being dynamically probed on the potential energy surface. 

At low energies a/k is small as expected, typical values being ~ 0.2 — 0.4. At highest 
energies reached, it is between 0..8 and 1.2 for various clusters with somewhat higher values 
for LJ clusters and smaller N. In the bulk system, a peak a/k k, 0.8 is reached at the melting 
point (from the solid phase) and then it remains nearly constant. The correspondence of 
the maximal roughness with melting point is very suggestive. One can visualize destruction 
of the crystal lattice being driven by large scale roughness of the potential. 


IV. CONCLUSION 


In this paper we have examined the behaviour of the largest Lyapunov exponent A as a 
function of energy in hnite clusters of 6-13 rare-gas and metal atoms, and in bulk rare-gas 
solids. These systems undergo a phase transition from a regime wherein the dynamics is 
purely oscillatory (involving individual particle vibrations) to a regime where the dynamics 
is both oscillatory as well as diffusive. 

Diffusive dynamics is linked to the presence of delocalized unstable modes in the bulk 
4^ . In small clusters the onset of the diffusion does not appear to enhance the chaoticity: 


the observed value of the Lyapunov exponent is smaller compared to the value expected 
by a simple extrapolation of the exponent from the low energy regime,namely from the 
oscillatory dynamics or the “solid” phase. This suppression of chaos, which we attribute to 
the correlated hopping dynamics, is strongest for smallest clusters but is then progressively 
reduced. It is possible that for particular clusters, the enhancing and suppressing effects of 
the unstable modes can balance and A(e) curve is smooth across the melting (infact LJ 17 
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shows no signature of melting according to this measure This conjecture can be tested 

by studies of the larger clusters with calculations of the participation ratios of the unstable 
modes which will clarify their role in the chaoticity of a dynamical system. 

One may intuitively expect that unstable modes (h e. negative curvatures) cause 
the dynamics to be chaotic. As noted by Dellago and Posch in their study of melting in 


two dimensional systems the fraction of unstable modes, which is a rough measure 
of negative curvatures, has a similar dependence on the parameters as A.These unstable 
modes can however become important only when particles are capable of large scale motion. 
(In a related context, it has been seen |2^ that A falls when a liquid is cooled through its 
glass-transition temperature, namely as the unstable modes get localized 

Using the framework of a geometric theory of Hamiltonian chaos, we compute an esti¬ 
mate for the Lyapunov exponent from the curvature of the potential energy surface and its 
fluctuation. We have studied the variation of these quantities with the temperature of the 
system and found that the mean curvature is always a monotonic function of energy but the 
variance has a simple energy dependence only for smaller clusters. In the coexistence region 
of 13 particle clusters—these are the cases in which potential energy surface has a deep 
global minimum which is well separated from the next lowest structure —a is nonmonotonic. 
The LJ bulk system shows singular behaviour for a at melting, which may indicate some 
sort of topological change in the conhguration space |T^ . 


The resulting estimate A is generally larger than A in the solid phase. For bulk, the 
discrepancy is small, but for clusters, the agreement is only qualitative. Curvature fluctu¬ 
ations for clusters are correlated, and this effectively reduces the parametric instability in 
the dynamics. The spectral features of the curvature fluctuations such as bandwidth are 
not well accounted for by the geometric theory even for the bulk system. At higher tem¬ 
peratures A is lower than A which is attributed to the unstable modes (negative curvatures) 
which are ignored in the geometric theory. The contribution of the unstable modes towards 
chaoticity (obtained by subtracting A from A and therefore only approximate) is small in 
the solid phase but can be as large as one-third in the liquid phase. Since a and k do not 
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show any singularity at melting for clusters, the parametric contribution coming from the 
change in topography of the potential energy surface is changing smoothly. Therefore, the 
fractional chaoticity coming from the unstable modes, 6Xu, also seems to vary continuously 
with energy. 

In summary, our application of the geometric theory to the dynamics of the melting 
transition for cluster and bulk systems has provided a satisfactory qualitative understanding 
of the underlying mechanisms in terms of the change in roughness of the potential energy 
surface, curvature fluctuations and and parametric instability. While agreement between 
theory and simulation is reasonable for the bulk system, for the case of hnite clusters, the 
situation is less satisfactory. The main source of the discrepancy seems to lie in the fact 
that in cluster systems, correlations are temporally long lived. This aspect requires to be 
incorporated within the present framework of the geometric theory (see e. g. Refs. [^p7|) 
in order to achieve quantitative accuracy. 
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FIGURES 

FIG. 1. (a) A as the function of reduced energy per particle for LJ clusters with ^^=6,7,9,11,13- 
Note that the list includes both magic and nonmagic clusters, (b) A for Gu clusters with Ai=6,7,13. 
The reduced energy scale is set by Ehuik (Eq. 2) 

FIG. 2. (a) Lyapunov exponents for the bulk LJ system of 256 particles in a cubical box with 
periodic boundary conditions. Circles refers to A and squares to the estimate A generated using 
Eq. 8 with r dehned in Eq. 10. Pluses (+) denote values of A calculated using T 2 = (b) 

Mean curvature (k), and fluctuation a for bulk LJ system as a function of energy, k and a are 
measured in units of frequency squared. 

FIG. 3. (a) Mean curvature k and (b) fluctuation a for LJat clusters with N=6, 11, 13 as a 
function of energy. Units are as in Fig. 2. 

FIG. 4. The estimate A for LJjv clusters with N=6, 11, 13 as the function of energy. Shown 
with circles are A calculated with r defined in Eq. 10; + are values of A calculated using T 2 = k^^'^ja. 
Also, for comparision, are the correspnding values of A from Fig 1 (squares) 

FIG. 5. (a) Mean curvature (fc), and (b) fluctuation (a) Guat clusters with N=6, 7, 13 as a 
function of energy, k and a are scaled by Eiyuik/mr^. 

FIG. 6 . The estimate A for Guat clusters with N=6, 7, 13 as a function of energy (same 
conventions as Fig. 4). 

FIG. 7. Power spectra of curvature distribution calculated along a trajectory for (a) bulk LJ 
(solid), (b) bulk LJ (liquid), (c) LJ 13 at energy= -2.4 corresponding to the temperature 34K 
(coexistence region) and, (d) Guia at e = —0.67 (liquid phase). The vertical axis is in arbitrary 
units. 
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TABLE 1. Typical timescales (in reduced units) associated with power spectra for various 


systems studied here. 


System 




' c 

A 

Gue solid 

0.30 

1.5 

0.8 

2.0 

0.2 

Liquid 

0.18 

0.4 

0.7 

1.1 

0.7 

Gui 3 solid 

0.28 

0.9 

0.8 

1.4 

0.3 

Liquid 

0.18 

0.36 

0.7 

0.9 

0.8 

LJe solid 

0.55 

1.6 

2.5 

10. 

0.05 

Liquid 

0.4 

0.8 

1.0 

1.0 

.25 

LJi 3 solid 

0.5 

1.4 

1.7 

2.5 

0.1 

Liquid 

0.3 

0.6 

1.0 

1.0 

0.3 

Bulk solid 

0.2 

0.5 

0.6 

0.6 

0.3 

Liquid 

0.2 

0.6 

0.4 

0.4 

0.8 


“ calculated with Eq. 10. 

^ calculated with Eq. 11 

Obtained from approximate upper cutoff of the power spectrum. 
^ Inverse bandwidth of the power spectrum (approximate) 
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